package fitness;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Random;

import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.fitting.GaussianFitter;
import org.apache.commons.math3.fitting.PolynomialFitter;
import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
import org.math.plot.utils.Array;

import bsh.EvalError;
import bsh.Interpreter;
import optimizers.commons.GAConstants;
import optimizers.ganeat.amplify.ConstantsAmply;
import optimizers.ganeat.math.ConstantsMath;
import util.MyEventHandlerAmplify;
import util.MyStepHandlerAmplify;
import util.MyUtil;
import model.Constants;
import model.Input;
import model.OligoSystemComplex;
import model.OligoSystemGeneral;

public class MathFitness extends Fitness implements Serializable {

	/**
         * 
         */
	private static final long serialVersionUID = 1L;
	private boolean debug = false;
	public static final double minvalue = 1e-2 * ConstantsAmply.amplitudeInput;

	public double findMax(double[] output) {
		double local = output[ConstantsAmply.maxDowntime]; // we might want to
															// let the system
															// settle for a
															// while, though
		for (int i = 1; i < output.length; i++) {
			if (output[i] > local) {
				local = output[i];
			} else if (output[i] < 0) {
				return -1;
			}
		}
		return local;
	}

	@Override
	public FitnessResult evaluateFitness(OligoSystemGeneral oligoModel) {
		try {
			MathResult result = new MathResult();
			HashMap<String, Double> errors = new HashMap<String, Double>();
			HashMap<String, Double> targets = new HashMap<String, Double>();
			HashMap<String, Double> actuals = new HashMap<String, Double>();
			boolean[] sequenceProtected = new boolean[oligoModel.nSimpleSequences];
			for (int i = 0; i < ConstantsMath.nInput; i++) {
				sequenceProtected[i] = true;
			}
			for (int i = ConstantsMath.nInput; i < oligoModel.nSimpleSequences; i++) {
				sequenceProtected[i] = false;
			}
			oligoModel.sequenceProtected = sequenceProtected;
			double[] reporterConcentration = new double[oligoModel.nSimpleSequences];
			for (int i = 0; i < oligoModel.nSimpleSequences; i++) {
				reporterConcentration[i] = 0;
			}
			reporterConcentration[ConstantsMath.nInput] = ConstantsMath.reporterConcentration;
			oligoModel.reporterConcentration = reporterConcentration;

			ArrayList<double[]> tests = ConstantsMath.logScale ? MyUtil
					.getTestsForMathLogScale(ConstantsMath.minTest,
							ConstantsMath.maxTest, ConstantsMath.nTests)
					: MyUtil.getTestsForMath(ConstantsMath.minTest,
							ConstantsMath.maxTest, ConstantsMath.nTests);

			double[][][] datas = oligoModel
					.getTimeSeriesComplexModelTillStableWithDifferentInputSet(tests);

			boolean allBelowThreshold = true;
			for (int j = 0; j < tests.size(); j++) {
				double[] inputs = tests.get(j);
				double actual;
				double[][] data = datas[j];
				int output = data.length - 1;
				if (output >= 0) {
					int firstStableTime = data[output].length - 1;
					if (firstStableTime == -1) {
						actual = -1;
					} else {
						while (data[output][firstStableTime] == 0
								&& firstStableTime > 0) {
							firstStableTime--;
						}
						if (firstStableTime >= Constants.maxTime - 1) {
							actual = -1;
						} else {
							actual = data[output][firstStableTime];
						}
					}
				} else {
					actual = -1;
				}
				if (actual < 0) {
					return minFitness(oligoModel);
				}
				if (actual > 1) {
					allBelowThreshold = false;
				}
				actuals.put(Array.toString(inputs), actual);
			}
			if (allBelowThreshold) {
				return minFitness(oligoModel);
			}

			if (!ConstantsMath.expression.startsWith("fit ")) {
				for (int j = 0; j < tests.size(); j++) {
					double[] inputs = tests.get(j);
					double target;
					StringBuilder builder = new StringBuilder();
					builder.append(ConstantsMath.preExpression);
					String expression = ConstantsMath.expression;
					for (int i = 0; i < ConstantsMath.nInput; i++) {
						expression = expression.replace("input" + i,
								String.valueOf(inputs[i]));
					}
					builder.append(expression);
					Interpreter interpreter = new Interpreter();
					try {
						interpreter.eval(builder.toString());
						Double r = (Double) interpreter.get("result");
						target = r;
					} catch (EvalError e) {
						e.printStackTrace();
						target = -1;
					}
					targets.put(Array.toString(inputs), target);
				}
			} else {
				String[] tmp = ConstantsMath.expression.split(" ");
				double[] targetFittedParams = new double[tmp.length - 2];
				double[] actualFittedParams = new double[tmp.length - 2];
				double[] best = null;
				switch (tmp[1]) {
				case "gaussian":
					GaussianFitter gaussianFitter = new GaussianFitter(
							new LevenbergMarquardtOptimizer());
					for (int j = 0; j < tests.size(); j++) {
						double[] inputs = tests.get(j);
						if (ConstantsMath.logScale) {
							gaussianFitter.addObservedPoint(
									Math.log(inputs[0]),
									actuals.get(Array.toString(inputs)));
						} else {
							gaussianFitter.addObservedPoint(inputs[0],
									actuals.get(Array.toString(inputs)));
						}
					}
					best = gaussianFitter.fit();
					for (int i = 0; i < tmp.length - 2; i++) {
						targetFittedParams[i] = Double.parseDouble(tmp[i + 2]);
						actualFittedParams[i] = best[i];
					}
					for (int j = 0; j < tests.size(); j++) {
						double[] inputs = tests.get(j);
						double target = 0;
						if (ConstantsMath.logScale) {
							target = actualFittedParams[0]
									* Math.exp(-(Math.log(inputs[0]) - actualFittedParams[1])
											* (Math.log(inputs[0]) - actualFittedParams[1])
											/ (2 * actualFittedParams[2] * actualFittedParams[2]));
						} else {
							target = actualFittedParams[0]
									* Math.exp(-(inputs[0] - actualFittedParams[1])
											* (inputs[0] - actualFittedParams[1])
											/ (2 * actualFittedParams[2] * actualFittedParams[2]));
						}
						targets.put(Array.toString(inputs), target);
					}
					break;
				case "square":
					PolynomialFitter fitter = new PolynomialFitter(
							new LevenbergMarquardtOptimizer());
					for (int j = 0; j < tests.size(); j++) {
						double[] inputs = tests.get(j);
						fitter.addObservedPoint(inputs[0],
								actuals.get(Array.toString(inputs)));
					}
					final double[] initSquare = { 0, 0, 1 };

					// Compute optimal coefficients.
					best = fitter.fit(new MyParametric(), initSquare);
					targetFittedParams[0] = Double.parseDouble(tmp[2]);
					actualFittedParams[0] = best[2];
					for (int j = 0; j < tests.size(); j++) {
						double[] inputs = tests.get(j);
						double target = actualFittedParams[0] * inputs[0]
								* inputs[0];
						targets.put(Array.toString(inputs), target);
					}
					break;
				case "plateau":
					for (int j = 0; j < tmp.length - 2; j++) {
						targetFittedParams[j] = Double.parseDouble(tmp[j + 2]);
					}
					final double[] initPlateau = { 1 };
					double bestMSE = Double.MAX_VALUE;
					for (int j = 0; j < tests.size() - 2; j++) {
						for (int k = j + 1; k < tests.size() - 1; k++) {
							fitter = new PolynomialFitter(
									new LevenbergMarquardtOptimizer());
							double[] arr1 = new double[tests.size()];
							double[] arr2 = new double[tests.size()];
							HashMap<String, Double> tmpTargets = new HashMap<String, Double>();
							double[] tmpActualFittedParams = new double[tmp.length - 2];
							for (int l = 0; l <= j; l++) {
								double[] inputs = tests.get(l);
								arr1[l] = actuals.get(Array.toString(inputs));
								arr2[l] = 0;
								tmpTargets.put(Array.toString(inputs), 0.0);
							}

							for (int l = j + 1; l <= k; l++) {
								double[] inputs = tests.get(l);
								arr1[l] = actuals.get(Array.toString(inputs));
								fitter.addObservedPoint(inputs[0], arr1[l]);
							}
							best = fitter
									.fit(new LineParametric(), initPlateau);
							double a = best[0];
							for (int l = j + 1; l <= k; l++) {
								double[] inputs = tests.get(l);
								arr2[l] = a;
								tmpTargets.put(Array.toString(inputs), a);
							}

							for (int l = k + 1; l < tests.size(); l++) {
								double[] inputs = tests.get(l);
								arr1[l] = actuals.get(Array.toString(inputs));
								arr2[l] = 0;
								tmpTargets.put(Array.toString(inputs), 0.0);
							}
							double mse = MyUtil.MSE(arr1, arr2);
							tmpActualFittedParams[0] = a;
							tmpActualFittedParams[1] = j;
							tmpActualFittedParams[2] = k - j;
							if (mse < bestMSE) {
								bestMSE = mse;
								targets = tmpTargets;
								actualFittedParams = tmpActualFittedParams;
							}
						}
					}
					break;
				}
				result.actualFittedParams = actualFittedParams;
				result.targetFittedParams = targetFittedParams;
			}

			double max_target = Double.MIN_VALUE;
			double max_actual = Double.MIN_VALUE;
			for (double target : targets.values()) {
				if (max_target < target) {
					max_target = target;
				}
			}
			for (double actual : actuals.values()) {
				if (max_actual < actual) {
					max_actual = actual;
				}
			}
			result.target_max = max_target;
			double scale = max_target / max_actual;
			if (scale > ConstantsMath.max_scale)
				scale = ConstantsMath.max_scale;
			if (scale < 1 / ConstantsMath.max_scale)
				scale = 1 / ConstantsMath.max_scale;

			for (String key : actuals.keySet()) {
				double error;
				double actual = actuals.get(key);
				double target = targets.get(key);
				if (actual == -1 || target == -1) {
					error = Double.MAX_VALUE;
				} else {
					if (!ConstantsMath.y_scaling) {
						error = (actual - target);
					} else {
						error = actual * scale - target;
					}
				}
				errors.put(key, error);
			}
			result.actuals = actuals;
			result.targets = targets;
			result.errors = errors;
			if (ConstantsMath.y_scaling) {
				result.scale = scale;
			}
			return result;
		} catch (Exception ex) {
			ex.printStackTrace();
			return minFitness(oligoModel);
		}
	}

	@Override
	public FitnessResult evaluateFitnessSimple(OligoSystemGeneral indiv) {
		// TODO Auto-generated method stub
		return null;
	}

	@Override
	public FitnessResult minFitness(OligoSystemGeneral indiv) {
		MathResult result = new MathResult();
		result.minFitness = true;
		return result;
	}
}

class MyParametric extends PolynomialFunction.Parametric {
	public double[] gradient(double arg0, double... arg1) {
		return new double[] { 0, 0, 1 };
	}
}

class LineParametric extends PolynomialFunction.Parametric {
	public double[] gradient(double arg0, double... arg1) {
		return new double[] { 1 };
	}
}